#include <iostream>
#include <fstream>
#include <vector>
#include <string>

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include<iostream>
#include<fstream>
#include<typeinfo>

#include "util.h"
#include "diffuse.h"

// 定义方阵的维度
#define nMAX 10  

using namespace std;


int main()
{

	char str1[64], str2[64], str3[64];

	sprintf(str1, "output_data");
	if (!access(str1, 0) == 0) mkdir(str1, 0777); // if  no "data" folder, construct it for output date

	char   buf_data[40];

	int ndim = 3;
 	
	if (ndim > nMAX) { cout << "矩阵阶数超过最大数nMAX !!!" << endl;  exit(-1); }

	// double  a[nMAX][nMAX] = { 0 };//定义一个nMAX*nMAX*的矩阵，用于存放数据
	double  **a = dmatrix(1, ndim, 1, ndim); // 或者我们可以动态开辟数组，这种方式记得在程序结束后释放内存。使用方式 free_dmatrix(a, 0, nMAX, 0, nMAX);
	
	// 怎样读取数据文件，比如data.txt中存储的是矩阵A数值
	ifstream infile;//定义读取文件流，相对于程序来说是in
	infile.open("data.txt");//打开文件
	for (int i = 1; i <= ndim; i++)//定义行循环
	{
		for (int j = 0; j <= ndim; j++)//定义列循环
		{
			infile >> a[i][j];//读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
		}
	}
	infile.close();//读取完成之后关闭文件


	ifstream infile2;
	infile2.open("data_b.txt");		// 打开文件
	for (int i = 1; i <= nMAX; i++) // 定义行循环
	{
		infile2 >> b[i]; // 读取一个值（空格、制表符、换行隔开）就写入到矩阵中，行列不断循环进行
	}
	infile2.close(); // 读取完成之后关闭文件


//////////////////////////////////// test  /////////////////////

	// Givens vector
	double x1 = 1;
	double x2 = 1;
	
	double *c, *s;

	Givens(x1, x2, &c, &s);


	printf("c = %f, s = %f \n ", *c, *s);

    return 0;
}

void Givens(double x1, double x2, double *c, double *s)
{
	// Given x=[x1, x2]^T in R^2, compute c and s such that Gx = [r, 0]^T, where r = |x|_2
	double tau;

	if (x2 == 0)
	{
		*c = sign(x1);
		*s = 0.0;
	}
	else
	{
		if (fabs(x2) > fabs(x1))
		{
			tau = x1 / x2;
			*s = sign(x2) / sqrt(1.0 + tau * tau);
			*c = (*s) * tau;
		}
		else 
		{
			tau = x2 / x1;
			*s = sign(x1) / sqrt(1.0 + tau * tau);
			*c = (*s) * tau;
		}
	}
	return;	
}


double sign(double x)
{
	if (x >= 0.0) 
		return 1.0; 
	else 
		return -1.0;
}


void print_vector(double *x, int n, char* buf)
{
	char buffer_x[32];
	
	FILE *fx;

	sprintf(buffer_x, "./%s/x.m", buf);

	fx = fopen(buffer_x, "w");
	for (int i = 1; i <= n; i++)
	{
		fprintf(fx, " %16.14f ", x[i]);
		fprintf(fx, "\n");
	}
	fclose(fx);

	return;
}



void print_matrix(double** phi, int n, char* buf)
{

	char buffer_phi[32];

	FILE *fphi;

	sprintf(buffer_phi, "./%s/a.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	return;
}


void print_matrix2(double** phi, double** eta, int n, char* buf)
{

	char buffer_phi[32], buffer_eta[32];

	FILE *fphi;
	FILE *feta;


	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);

	return;
}



void print_matrix3(double** phi, double** eta, double **xx,  int n,  char* buf)
{

	char buffer_phi[32], buffer_eta[32], buffer_xx[32];

	FILE *fphi;
	FILE *feta;
	FILE *fxx;

	sprintf(buffer_phi, "./%s/a.m", buf);
	sprintf(buffer_eta, "./%s/L.m", buf);
	sprintf(buffer_xx,  "./%s/U.m", buf);

	fphi = fopen(buffer_phi, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fphi, " %16.14f ", phi[i][j]);
		}
		fprintf(fphi, "\n");
	}
	fclose(fphi);

	feta = fopen(buffer_eta, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(feta, " %16.14f ", eta[i][j]);
		}
		fprintf(feta, "\n");
	}
	fclose(feta);


	fxx = fopen(buffer_xx, "w");
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			fprintf(fxx, " %16.14f ", xx[i][j]);
		}
		fprintf(fxx, "\n");
	}
	fclose(fxx);

	return;
}


